* VAR code for results in Leduc, S. and K. Sill, "Expectations and Economic Fluctuations: An Analysis Using Survey Data,"
* Review of Economic and Statistics 95(4), October 2013, 1352-1367.
* This code generates the responses to a positive shock to expected unemployment in the first column of Figure 3, which uses the Livingston expectation data.
* To reproduce the results in Figure 3, the positive expectation shock must be converted to a negative one.

********************************************************************************
*                                                                              *
* Note:  Computes 3 Confidence Intervals Simultaneously (90%, 80%, 68%)        *
*                                                                              *
********************************************************************************

********************************************************************************
* Note:  This program renormalizes the impulse responses to VAR shocks
*        by dividing all impulse responses by the own-contemporaneous response
*        (see the two "Patch Sections")
*
*        This program also computes historical decompositions & structural shock
*        estimates.
********************************************************************************

*********************************************************************************************
* Program Purpose:
* (1) Reads in semiannual Livingston survey data:
*     yy:1 = the June survey (which forecasts DEC)
*     yy:2 = the Dec  survey (which forecasts JUNE)
*
* (2) Reads in monthly macro variables as semiannual data and constructs all transformations:
*     yy:1 = avg of May(5,t), June(6,t), July(1,t+1), Aug(2,t+1), Sep(3,t+1), and Oct(4,t+1)
*     yy:2 = avg of Nov(5,t), Dec(6,t),  Jan(1,t+1),  Feb(2,t+1), Mar(3,t+1), and Apr(4,t+1)
*
* (3) VAR: pdot, expec, u, r
*     IDENTIFICATION:  Cholesky schemes
*
* (4) Computes optimal lag selection via AIC, SIC, but also allows for arbitrary lag choice
*
* (5) Estimates VAR
*
* (6) Computes impulse responses, bootstrap confidence intervals for VAR shocks, with all
*     impulse responses renormalized by dividing by the contemporaneous own response
*
*********************************************************************************************

cal 1947 1 2
alloc 0 2011:3

****************************************
*                                      *
* Part 1. Read Livingston Survey Data  *
*                                      *
****************************************

open data c:\myfiles\papers\news\code\data_restat.rat
data(format=rats,org=obs) 1947:1 2011:3 ruc6med

****************************************************************
*                                                              *
* Part 2. Input Macro Variables                                *
*         As Semi-Annual Observations                          *
*                                                              *
****************************************************************


*******************************
* 1. Jan & July Monthly Values
*******************************
data(format=rats,org=obs,select=1) 1947:1 2011:3  ruc  rtb3
set ruc1      / = ruc(t)
set r1        / = rtb3(t)

clear ruc rtb3

*******************************
* 2. Feb & Aug Monthly Values
*******************************
data(format=rats,org=obs,select=2) 1947:1 2011:3  ruc  rtb3
set ruc2      / = ruc(t)
set r2        / = rtb3(t)

clear ruc rtb3

*******************************
* 3. Mar & Sept Monthly Values
*******************************
data(format=rats,org=obs,select=3) 1947:1 2011:3  ruc  rtb3
set ruc3      / = ruc(t)
set r3        / = rtb3(t)

clear ruc rtb3

*******************************
* 4. Apr & Oct Monthly Values
*******************************
data(format=rats,org=obs,select=4) 1947:1 2011:3  ruc  rtb3 cpi
set ruc4      / = ruc(t)
set r4        / = rtb3(t)
set p4        / = cpi(t)

clear ruc rtb3 cpi

*******************************
* 5. May & Nov Monthly Values
*******************************
data(format=rats,org=obs,select=5) 1947:1 2011:3  ruc  rtb3
set ruc5      / = ruc(t)
set r5        / = rtb3(t)

clear ruc rtb3

*******************************
* 6. Jun & Dec Monthly Values
*******************************
*data(format=rats,org=obs,select=6) 1948:1 2010:12  ruc  rtb3
data(format=rats,org=obs,select=6) 1947:1 2011:3  ruc  rtb3
set ruc6      / = ruc(t)
set r6        / = rtb3(t)

clear ruc rtb3


*************************************************************************************
*                                                                                   *
* Part 3. Compute VAR Semi-Annual Macro Data                                        *
*                                                                                   *
* Note 1. The first period  in each year (yyyy:1) is the Livingston June Survey     *
*         The second period "   "     "  (yyyy:2)  "  "     "       Dec  Survey     *
*                                                                                   *
*************************************************************************************

print(dates) / ruc5
* Unemployment Rates
set urate   /  = (ruc5(t)  + ruc6(t)  + ruc1(t+1)  + ruc2(t+1)  + ruc3(t+1)  + ruc4(t+1))/6.0
*print(dates) / urate

* Interest Rates
set r       /  = (r5(t)    + r6(t)    + r1(t+1)    + r2(t+1)    + r3(t+1)    + r4(t+1))/6.0

* Inflation Rates
set inf     /  = (log(p4(t+1)/p4(t)  ))*2.0*100.0


*Oil shock quantitative dummy
set oshock  /  = 0.0
set oshock 1973:2 1973:2 =  7.8
set oshock 1978:2 1978:2 =  8.9
set oshock 1980:1 1980:1 =  7.2
set oshock 1990:1 1990:1 =  8.8

set oshockreplace  / = oshock

*Fiscal policy dummy
set gshock  /  = 0.0
set gshock 1965:1 1965:1 = 1
set gshock 1980:1 1980:1 = 1
set gshock 2001:2 2001:2 = 1

clear  ruc1     ruc2     ruc3    ruc4    ruc5    ruc6
clear  r1       r2       r3      r4      r5      r6
clear  p4


************************************************************
*                                                          *
* Part 4.  Program Parameters For Lag Selection Routine    *
*                                                          *
************************************************************

*******************************************
* Set Generic Variables to Use in the VAR
*******************************************

set pdot      /  = inf;  		* actual inflation measure
set expec     /  = ruc6med; 		* expected unemployment rate
set u         /  = urate; 		* unemployment rate measure
set r         /  = r;  			* nominal rate measure

*print(dates) / pdot expec u r


*********************************************************************************************
* Set Sample Periods & Other Parameters
*
* Note 1. We do not want lags of RHS variables to extend into the period before 1952:1
*         (in the first sample) or into the the period before 1979:2 (in the second sample).
*         Thus, given the samples listed below, we adjust the starting point of each sample
*         in the estimation commands, depending on the number of lags used in the VAR.
*********************************************************************************************

compute sbeg1 = 1961:1      ;    * first period for estimation
compute send1 = 2008:2      ;    * last period  for estimation

declare rectangular[integer]  lagchoice(1,3)

        * Column 1  :   Corresponds to AIC (determined below)
        * Column 2  :   Corresponds to SIC (determined below)
        * Column 3  :   Any number you want

compute lagchoice(1,3) =  2 ; * enter an integer for number of lags to use, for use when you do
                              * not want to use the AIC or SIC (pre-Volcker period)

compute LS1   = 3;          ; * Lag Scheme(1 = AIC, 2 = SIC, 3 = other)

compute nvar  = 4           ; * number of variables in the VAR

                              * List order of variables as given below in the VARIABLES command

compute [vector[string]] varnames = ||'pdot', 'expec', 'u', 'r'||


******************************************************
*                                                    *
* Part 5.  Compute AIC & SIC Minimizing Lag Lengths  *
*                                                    *
******************************************************

set oshock1 / = oshock(t-1)
set oshock2 / = oshock(t-2)

set gshock1 / = gshock(t-1)
set gshock2 / = gshock(t-2)


source(noecho) c:\myfiles\papers\news\code\aicsic.src


@aicsic(more, nosuppress) sbeg1  send1  4  lagchoice(1,1) lagchoice(1,2)
# pdot expec u r
*# constant oshock oshock1 oshock2
# constant


*display;display;


                 ***************************************
                 *                                     *
                 *         Begin Bootstrap Routine     *
                 *                                     *
                 ***************************************

*********************************************************************************************
* This program produces the following:
*   IRPOINT_ = NVARS x NSHOCKS RECT[SERIES] of impulse response point estimates,
*             where the response of the ith variable to a shock to the jth equation is
*             IRPOINT(i,j)
*   LBAND_, UBAND_ =
*             NVARS x NSHOCKS RECT[SERIES] of lower and upper points of the confidence
*             interval, each dimensioned in the same way as IRPOINT
* Thus, the impulse response and confidence intervals for y(i) in response to a shock to
* equation j are given by IRPOINT_(i,j) and LBAND_(i,j), UBAND_(i,j).
*
* The (i,j) refer to the order in which the variables are listed in Part 4.a - 4.b (which is
* the same order that you must construct the VAR in Part 6), not the Cholesky ordering defined
* in Part 5.
*
*
*
* Note:  all bootstrap parameters and data types use "_" as the last character of the name
*
*
**********************************************************************************************

**********************************
*                                *
* Bootstap Parameters            *
*                                *
**********************************

compute sbeg1_     = sbeg1+lagchoice(1,LS1); * first period to use to estimate the VAR
                                             * (all lagged values MUST be defined at this date)

compute send1_     = send1;                  * last  period to use to estimate the VAR

compute subtitle_  = %datelabel(sbeg1_)+'-'+%datelabel(send1_)

compute ndraws_    = 1000;   * number of draws

compute nvars_     = nvar+1; * number of variables in the model = VAR eqns + identities
                             * (you must include at least one identity, thus nvars > nshocks)

compute nshocks_   = nvar;      * number of shocks in the VAR (i.e., nvars minus number of identities)
                             * (this is the number of equations in the VAR, not incl. identities)

compute nlags_     = lagchoice(1,LS1);   * number of lags in the VAR

compute nirsteps_  = 25;     * number of impulse response steps to compute

compute width1_     = 0.90;   * enter width for confidence interval 1 (0.95 means 95% intervals)
compute width2_     = 0.80;   * enter width for confidence interval 2 (0.95 means 95% intervals)
compute width3_     = 0.68;   * enter width for confidence interval 3 (0.95 means 95% intervals)


compute [vector] percentiles1_ $
                    = ||(1.0-width1_)/2.0,(1.0+width1_)/2.0||;  * percentiles for boostrapped error bands 1

compute [vector] percentiles2_ $
                    = ||(1.0-width2_)/2.0,(1.0+width2_)/2.0||;  * percentiles for boostrapped error bands 2

compute [vector] percentiles3_ $
                    = ||(1.0-width3_)/2.0,(1.0+width3_)/2.0||;  * percentiles for boostrapped error bands 3


**********************************
*                                *
* Dimension Objects              *
*                                *
**********************************


declare vector[series]        path_(nshocks_) actser_(nvars_)
declare rectangular           pathmatrix_(send1_-sbeg1_+1,nshocks_) irall_(ndraws_,nirsteps_*nvars_) $
                                                                    iroil_(ndraws_,nirsteps_*nvars_)

declare rectangular[series]   ssim_(ndraws_,nvars_) lband1_(nvars_,nshocks_) uband1_(nvars_,nshocks_)  $
                                                    lband2_(nvars_,nshocks_) uband2_(nvars_,nshocks_)  $
                                                    lband3_(nvars_,nshocks_) uband3_(nvars_,nshocks_)  $
                                                    irpointoil_(nvars_,1)                              $
                                                    lbandoil1_(nvars_,1)     ubandoil1_(nvars_,1)      $
                                                    lbandoil2_(nvars_,1)     ubandoil2_(nvars_,1)      $
                                                    lbandoil3_(nvars_,1)     ubandoil3_(nvars_,1)



declare vector[integer]       sbegs_(nvars_)




*****************************************************************
*                                                               *
* Place All Relevant Data Into VEC[SER] = ACTSER_               *
*                                                               *
*   Note 1. The ordering does not affect the Cholesky Decomp    *
*   Note 2. The ordering only affects the position of the VAR   *
*           equations in the model                              *
*   Note 3. You must have nvars entries, VAR equation variables *
*           first, then variables defined by identity           *
*                                                               *
*****************************************************************

* 4.A.  VAR Equation Variables (these are the dependent variables of the VAR)

set actser_(1)    /  = pdot
set actser_(2)    /  = expec
set actser_(3)    /  = u
set actser_(4)    /  = r


do i = 1, nshocks_
   labels actser_(i)
   # varnames(i)
end do i


* 4.B.  Variables Defined By Identity (these are the dependent variables in the identities)

set realrate      /  = r(t) - pdot(t)

set actser_(5)    /  = realrate(t)

* 4.C. Create Vector Of Starting Dates For Each ACTSER_(i)

do i = 1,nvars_
   inquire(series=actser_(i)) sbegs_(i) scrap
end do i


***************************************
*                                     *
* Define Cholesky Ordering            *
*                                     *
***************************************

compute [vector[integer]] corder_ = ||2,3,1,4||
compute choleskystring_ = ' / Order = expec, Pdot, U, R'

**********************************************************
*                                                        *
* Define The VAR                                         *
*                                                        *
*    Note 1. List variables in the order given in 4.A    *
*                                                        *
**********************************************************


system(model=VAR1eqns_)
     variables  actser_(1) actser_(2) actser_(3) actser_(4)
     lags 1 to nlags_
     det constant oshock{0 to nlags_} gshock{0 to nlags_}
end(system)



***********************************************
*                                             *
* Define All Identities                       *
*                                             *
*                                             *
* Note: you must hace at least one identity   *
***********************************************



* Real Rate

equation(identity) realrate_id   actser_(5)
# actser_(4)  actser_(1)
associate          realrate_id
# 1.0        -1.0



*********************************************************************
*                                                                   *
* Estimate The VAR Model                                            *
*         Capture Residuals                                         *
*         Capture Coefficients                                      *
*         Compute Cholesky Factor                                   *
*         Group Estimated VAR Equations & Identities Into a Model   *
*                                                                   *
*********************************************************************

estimate(residuals=resid_,coeffs=beta_)  sbeg1_  send1_

*linreg expec sbeg1_ send1_ resids
*# constant expec{1 to 2}
*statistics resids


*linreg expec sbeg1_ send1_ resids
*# constant expec{1 to 2} pdot{1 to 2} u{1 to 2} r{1 to 2} oshock{0 to 2} gshock{0 to 2}
*# constant expec{1 to 2} pdot{1 to 2} u{1 to 2} r{1 to 2}
*statistics resids
*exclude(title="Granger Causality Test")
*# pdot{1 to 2} u{1 to 2} r{1 to 2} oshock{0 to 2} gshock{0 to 2}
**# pdot{1 to 2} u{1 to 2} r{1 to 2}

compute cfactor_ = %psdfactor(%sigma,corder_)

group ID1model_  realrate_id

compute [model] model1_     = VAR1eqns_ + ID1model_


*******************************************************************************
*                                                                             *
* Compute & Capture Impulse Response Point Estimates, Variance Decomps,       *
* Historical Decompositions & Structural Shocks                               *
*                                                                             *
* Note 1.  The impulse response point estimates are in the NVARS_ x NSHOCKS_  *
*          RECT[SERIES] called IRPOINT_(i,j), i = 1,...,NVARS_ and            *
*          j = 1,...,NSHOCKS_.  The response of variable i to a shock         *
*          to equation j is IRPOINT_(i,j), where the number is that given in  *
*          Part 4A-4B, not by the Cholesky chosen in Part 5                   *
*                                                                             *
*******************************************************************************

display 'Cholesky Factor = ' cfactor_

****************************************
* Compute Structural Errors
*    = S_ERROR, a nshocks_ x 1 VEC[SER]
****************************************

declare vector[integer]  residlist
declare vector[series]   s_error(nshocks_)

enter(varying) residlist
# resid_(1)
do i = 2, nshocks_
   enter(varying) residlist
   # residlist resid_(i)
end do i

make rf_errormat  sbeg1_  send1_
# residlist

compute s_errormat = (rf_errormat)*(tr(inv(cfactor_)))

do i = 1, nshocks_
   set s_error(i) sbeg1_ send1_ = s_errormat(t-(sbeg1_-1),i)
end do i



*******************************************************************
* Compute: Var Decomp, Hist Decomp, Impulse Responses
* Note: HISTORY command will not work on a model with identities
*******************************************************************

errors(model =model1_,decomp=cfactor_)                      *  nirsteps_  *
history(model=VAR1eqns_,decomp=cfactor_,results=hist_,noadd)  *  (send1_-sbeg1_+1) sbeg1_
impulse(model=model1_,decomp=cfactor_,results=irpoint_)     *  nirsteps_  *

*print(dates) / hist_

******************************************************************************************
*         PATCH:  divide all impulse responses by contemporaneous own response
******************************************************************************************
do j = 1, nshocks_
   compute rescale = irpoint_(j,j)(1)
   do i = 1, nvars_
      set irpoint_(i,j) 1 nirsteps_ = irpoint_(i,j)/rescale
   end do i
end do j
******************************************************************************************



*******************************************************************************
*                                                                             *
* Compute & Capture Impulse Response Point Estimates to Oil Dummy Shock       *
*                                                                             *
* Note:  we compute this as the difference between 2 forecasts, one of which  *
*        has a (temporary) unit increase in Hamilton's oil dummy measure in   *
*        the first period of the forecast.                                    *
*                                                                             *
*******************************************************************************

****************************
* Compute Baseline Forecast
****************************

forecast(model=model1_,results = noshockfor_) * nirsteps_  sbeg1_

***********************************************
* Compute Forecast With Higher Dummy Variable
***********************************************

set oshock sbeg1_ sbeg1_ = oshock + 1.0

forecast(model=model1_,results = shockfor_)   * nirsteps_  sbeg1_

set oshock sbeg1_ sbeg1_ = oshockreplace


*************************************
* Compute Implied Impulse Responses (as the difference between the forecasts)
*************************************

do i = 1, nvars_

   set irpointoil_(i,1) / = shockfor_(i) - noshockfor_(i)

end do i

clear shockfor_  noshockfor_



*******************************************************************************
*                                                                             *
* Compute Simulated VAR Data                                                  *
*                                                                             *
* Note 1. The simulated values are in the NDRAWS x NVARS RECT[SERIES]         *
*         called SSIM_(draw,i),where the row (draw) indicates the replication *
*         number and the column (i) indicates the variable simulated          *
*                                                                             *
*******************************************************************************

* Conduct Simulation

smpl sbeg1_ send1_
do draw_ = 1, ndraws_

   boot entries_ / sbeg1_ send1_

   do i = 1,nshocks_
      set path_(i) / = resid_(i)(entries_(t))
   end do i

   do i = 1, nshocks_
      do j = 1, send1_-sbeg1_+1
         compute pathmatrix_(j,i) = path_(i)(sbeg1_-1+j)
      end do j
   end do i

   forecast(model=model1_,matrix=pathmatrix_, results=results_)

   do i = 1,nvars_
      set ssim_(draw_,i) / = results_(i)
   end do i

end do draw_


smpl
* Reset SMPL Above This Line


*********************************************************************************************
*                                                                                           *
* Estimate VAR  on Simulated Data                                                           *
*          Compute Impulse Responses                                                        *
*          Calculate CIs                                                                    *
*                                                                                           *
* Note 1.  The lower and upper confidence intervals are given in the NVARS x NSHOCKS        *
*          RECT[SERIES] called LBAND_(var,shock) and UBAND_(var,shock), where the row (var) *
*          indicates the response variable and the column (shock) indicates the equation    *
*          that is shocked. Note that the variable number and equation-shock number are     *
*          given in Part 4.A-4.B, not the Cholesky ordering in Part 5.                      *
*                                                                                           *
*********************************************************************************************

* Splice Initial Conditions (given by history) On Top Of Simulated Values

do draw_ = 1, ndraws_
   do var_ = 1, nvars_

      set ssim_(draw_,var_) sbegs_(var_) sbeg1_-1  = actser_(var_)

   end do var_
end do draw_


* Impulse Responses

do shock_ = 1,nshocks_

do draw_=1,ndraws_

   * Estimation

   system(model=VAR2eqns_)
     variables ssim_(draw_,1) ssim_(draw_,2) ssim_(draw_,3) ssim_(draw_,4)
     lags 1 to nlags_
     *det constant oshock{0 to nlags_}
     det constant oshock{0 to nlags_} gshock{0 to nlags_}
     *det constant
   end(system)

   estimate(noprint,noftests,coeffs=simbeta_) sbeg1_ send1_



   * Compute Coefficient Averages (over draws)[Not Used In This Program]

   if draw_.eq.1
     {
       compute musimbeta_ = (1.0/ndraws_)*simbeta_
     }
   else
     {
       compute musimbeta_ = musimbeta_ + (1.0/ndraws_)*simbeta_
     }

   * Define Identities

   * Real Rate
   equation(identity) realrate_id   ssim_(draw_,5)
   # ssim_(draw_,4)  ssim_(draw_,1)
   associate          realrate_id
   # 1.0           -1.0


   group ID2model_  realrate_id

   compute [model] model2_     = VAR2eqns_ + ID2model_


  * Compute Impulse Responses

   * Note 1: the equation shocked by the loop index, called "shock", is given by the order of
   *         equations in the VAR, not the order given by the cholesky factorization
   * Note 2: IR is an Nvars x 1 RECTANG[SER], with the ordering given by the ordering of the
   *         equations in the VAR, not the order of the cholesky factorization, and the
   *         ordering of the identities


   compute cfactor_ = %psdfactor(%sigma,corder_)
   impulse(model=model2_,noprint,decomp=cfactor_,results=ir_) * nirsteps_ shock_




   *****************************************************************************************
   *        PATCH:  divide all impulse responses by contemporaneous own response
   *****************************************************************************************
   compute rescale = ir_(shock_,1)(1)
   do i = 1, nvars_
     set ir_(i,1) 1 nirsteps_ = ir_(i,1)/rescale
   end do i
   *****************************************************************************************




   * Capture Impulse Responses For This Draw


   do var_ = 1,nvars_
      do step_ = 1,nirsteps_
          compute irall_(draw_,step_+(var_-1)*nirsteps_) = ir_(var_,1)(step_)
      end do step_
   end do var_

  * Compute & Capture Impulse Responses For Oil Dummy

   if shock_.eq.1
   {

      * Compute Baseline Forecast

      forecast(model=model2_,results = noshockfor_) * nirsteps_  sbeg1_

      * Compute Forecast With Higher Dummy Variable

      set oshock sbeg1_ sbeg1_ = oshock + 1.0

      forecast(model=model2_,results = shockfor_)   * nirsteps_  sbeg1_

      set oshock sbeg1_ sbeg1_ = oshockreplace

      * Compute Impulse Response As Difference Between The Forecasts & Capture

      do var_ = 1,nvars_
         do step_ = 1,nirsteps_
            compute iroil_(draw_,step_+(var_-1)*nirsteps_) = shockfor_(var_)(sbeg1_-1+step_)   $
                                                           - noshockfor_(var_)(sbeg1_-1+step_)
         end do step_
      end do var_

      clear noshockfor_ shockfor_

   }

end do draw_


   * Compute Confidence Intervals For This VAR Shock & The Oil Dummy Shock (by variable and step)

   do var_ = 1,nvars_
      do step_ = 1,nirsteps_

         * VAR Shock Responses
         compute col_ = %xcol(irall_,step_+(var_-1)*nirsteps_)
         compute lowerupper1_ = %fractiles(col_,percentiles1_)
         compute lowerupper2_ = %fractiles(col_,percentiles2_)
         compute lowerupper3_ = %fractiles(col_,percentiles3_)


         set lband1_(var_,shock_) step_ step_ = lowerupper1_(1)
         set lband2_(var_,shock_) step_ step_ = lowerupper2_(1)
         set lband3_(var_,shock_) step_ step_ = lowerupper3_(1)


         set uband1_(var_,shock_) step_ step_ = lowerupper1_(2)
         set uband2_(var_,shock_) step_ step_ = lowerupper2_(2)
         set uband3_(var_,shock_) step_ step_ = lowerupper3_(2)


         * Oil Dummy Shock Responses
         if shock_.eq.1
         {
            compute coloil_ = %xcol(iroil_,step_+(var_-1)*nirsteps_)

            compute lowerupperoil1_ = %fractiles(coloil_,percentiles1_)
            compute lowerupperoil2_ = %fractiles(coloil_,percentiles2_)
            compute lowerupperoil3_ = %fractiles(coloil_,percentiles3_)

            set lbandoil1_(var_,1) sbeg1_-1+step_   sbeg1_-1+step_  =  lowerupperoil1_(1)
            set lbandoil2_(var_,1) sbeg1_-1+step_   sbeg1_-1+step_  =  lowerupperoil2_(1)
            set lbandoil3_(var_,1) sbeg1_-1+step_   sbeg1_-1+step_  =  lowerupperoil3_(1)

            set ubandoil1_(var_,1) sbeg1_-1+step_   sbeg1_-1+step_  =  lowerupperoil1_(2)
            set ubandoil2_(var_,1) sbeg1_-1+step_   sbeg1_-1+step_  =  lowerupperoil2_(2)
            set ubandoil3_(var_,1) sbeg1_-1+step_   sbeg1_-1+step_  =  lowerupperoil3_(2)

         }

      end do step_
   end do var_

end do shock_



*********************************************
*                                           *
* Graphs of Impulse Responses:  VAR Shocks  *
*                                           *
*********************************************


*********************
* Expectations Shock
*********************

spgraph(vfields=3,hfields=2,header='Expectations Shock')

   graph(header='Inflation Response',nodates,number=0) 7
   # lband1_(1,2)     1 nirsteps_   1
   # irpoint_(1,2)    1 nirsteps_   4
   # uband1_(1,2)     1 nirsteps_   1
   # lband2_(1,2)     1 nirsteps_   3
   # uband2_(1,2)     1 nirsteps_   3
   # lband3_(1,2)     1 nirsteps_   6
   # uband3_(1,2)     1 nirsteps_   6

   graph(header='Expectations Response',nodates,number=0) 7
   # lband1_(2,2)     1 nirsteps_   1
   # irpoint_(2,2)    1 nirsteps_   4
   # uband1_(2,2)     1 nirsteps_   1
   # lband2_(2,2)     1 nirsteps_   3
   # uband2_(2,2)     1 nirsteps_   3
   # lband3_(2,2)     1 nirsteps_   6
   # uband3_(2,2)     1 nirsteps_   6

   graph(header='Unemployment Response',nodates,number=0) 7
   # lband1_(3,2)     1 nirsteps_   1
   # irpoint_(3,2)    1 nirsteps_   4
   # uband1_(3,2)     1 nirsteps_   1
   # lband2_(3,2)     1 nirsteps_   3
   # uband2_(3,2)     1 nirsteps_   3
   # lband3_(3,2)     1 nirsteps_   6
   # uband3_(3,2)     1 nirsteps_   6

   graph(header='Interest Rate Response',nodates,number=0) 7
   # lband1_(4,2)     1 nirsteps_   1
   # irpoint_(4,2)    1 nirsteps_   4
   # uband1_(4,2)     1 nirsteps_   1
   # lband2_(4,2)     1 nirsteps_   3
   # uband2_(4,2)     1 nirsteps_   3
   # lband3_(4,2)     1 nirsteps_   6
   # uband3_(4,2)     1 nirsteps_   6

spgraph(done)




















































































































































































































































































































